RMFinventory vignette

Introduction

This vignette presents the sampling and ABA modeling approaches that were used for the enhanced forest inventory (EFI) pilot project at Romeo Malette Forest (RMF) using Leica SPL100 lidar (SPL) data collected wall-to-wall during the summer of 2018. The package (and this vignette) are written for the RMNF project but the intent is to make the overall approach applicable to other forest management units.

First, we will perform a structurally-guided stratified sampling to select locations for new plots establishement. Then, we will detail how plot-level forest attribute data can be calculated or modeled from published allometric equations. Finally, we will use the field inventory data and the SPL point cloud to build random forest predicitve models and generate wall-to-wall maps of forest attributes.

We start by loading the packages that we are going to use

library(RMFinventory)
library(raster)
library(sf)
library(dplyr)
library(RStoolbox)
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 3.6.3

Select candidate cells from inventory layers

For this section we assume that wall-to-wall lidar metrics are available:

wall_metrics <- brick(system.file("extdata","wall_metrics_small.tif", package = "RMFinventory"))
names(wall_metrics) <- c("avg", "cov", "std","p10", "p20","p50","p70","p95", "p99","d0","d2","d4","dns")

#Coordinate system of wall-to-wal raster
wall_crs <- raster::crs(wall_metrics)

We will only select cells that are suitable for plot establishment using the followqign process:

1. Open a shapefile that contains the areas to keep/mask as a sf object
2. Perform some operations such as querying some attributes or create buffers
3. Convert the sf object as a sp object for cpompatibility with the raster package
4. USe the mask function of the raster package that will assign NA values to all cells not covered (i.e. centroid not inside polygon) by the sp object. ## PCI/BMI polygons

poly <- st_read(system.file("extdata/inventory_polygons","inventory_polygons.shp", package = "RMFinventory"))
#> Reading layer `inventory_polygons' from data source `C:\Users\queinnec.stu\Documents\R\win-library\3.6\RMFinventory\extdata\inventory_polygons\inventory_polygons.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 632 features and 111 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 431100 ymin: 5337700 xmax: 438560 ymax: 5343240
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=17 +ellps=GRS80 +units=m +no_defs

#Match CRS of wall-to-wall raster layer
poly <- st_transform(poly, wall_crs@projargs)
plot(poly["POLYTYPE"], axes = TRUE)

plot(poly["OWNER"], axes = TRUE)

poly_subset <- poly[poly$POLYTYPE == "FOR" & poly$OWNER == 1, ]
plot(st_geometry(poly_subset), axes = TRUE, col = "red")

poly_subset <- st_union(poly_subset)
plot(st_geometry(poly_subset), axes = TRUE, col = "red")

#Transform to SpatialPolygonDataFrame object for compatibility with the raster package
poly_subset <- sf::as_Spatial(poly_subset)

#Mask wall-to-wall layer
wall_poly <- raster::mask(wall_metrics, mask = poly_subset)
plot(wall_poly$avg)

Road layer

# Open roads layer, select suitable road types and create a buffer 

roads <- st_read(system.file("extdata/roads","roads.shp", package = "RMFinventory"))
#> Reading layer `roads' from data source `C:\Users\queinnec.stu\Documents\R\win-library\3.6\RMFinventory\extdata\roads\roads.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 167 features and 36 fields
#> geometry type:  MULTILINESTRING
#> dimension:      XY
#> bbox:           xmin: 431100 ymin: 5337700 xmax: 438560 ymax: 5343240
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=17 +ellps=GRS80 +units=m +no_defs
roads <- st_transform(roads, wall_crs@projargs)
plot(roads["RDTYPE"])

# Select only suitable road types
# Highway (H), Municipal (M), Primary (P), Branch (B), Clay/mineral surface haul (C) and roads accessed when dry or frozen (d) , graveled (g), industrial grade road (i), highway or municipal road(r), yearly accessible (y)

roads_subset <- roads[roads$RDTYPE %in% c("H", "M", "P", "B", "C") &
                        roads$RDMOD %in% c("d", "g", "i", "r", "y"), ]

# Dissolve roads layer to work with only 1 feature
roads_subset <- st_union(roads_subset)
plot(st_geometry(roads_subset))

#Create buffers of 30 m and 200 m
roads_buffer_30m <- st_buffer(roads_subset, dist = 30)
roads_buffer_200m <- st_buffer(roads_subset, dist = 200)

#Take the symetrical difference between buffers to keep only roads >- 30 m AND <= 200 m
roads_buffer <- st_sym_difference(roads_buffer_200m, roads_buffer_30m)
plot(st_geometry(roads_buffer))

#Transform to SpatialPolygonDataFrame object for compatibility with the raster package
roads_buffer <- sf::as_Spatial(roads_buffer)

#Mask wall-to-wall layer
wall_poly_roads <- raster::mask(wall_poly, mask = roads_buffer)
plot(wall_poly_roads$avg)

Additional masking layers

In this example we only mask by inventory polygons and distance to road but other layers can be considered (such as past or planned depletions). The same approach as the one described above can be used. # Stratified random sampling of cells for new plots establishement

When generating forest inventory attributes using an area-based-approach (ABA), the ground plots data used to calibrate ABA regression models need to be representative as much as possible of the full range of forest structure variability within the study area. If this is not the case, regression models might perform poorly in underrepresented forest types (White et al., 2013). Lidar metrics such as height percentiles, cover or height variability can be used to design a sampling network driven by forest structure.

Principal Component Analysis (PCA) is a method used to summarize the variability of a large number of highly correlated LiDAR structural metrics into a smaller number of uncorrelated variables. The feature space created by the generated principal components can then be stratified into classes that will represent specific types of forest structural conditions. Random sampling can then be performed within each of these classes to ensure a representative characterization of all forest structures occurring across the study area.

Calculate PCA values of candidate cells

We assume that a set of lidar metrics have been calculated across the entire study area and that candidate cells for new plots establishment have been determined using criteria such as distance and access to roads, forested land cover etc.

Note that the resolution of the lidar metrics maps is determined based on the size of the plots that will be established. It is important the the area of each cell (pixel) corresponds to the area of the plots. For fixed-radius circular plots of 11.28 m radius, corresponding to 400 m2, the resolution should be 20 m x 20 m.

In the example of RMF, we selected only cells located between 30 m and 200 m from roads and forest inventory polygons determined from photo interpretation, were used to select only cells considered to be forested.

# Plot wall-to-wall
plot(wall_metrics$p95, zlim = c(0, 30), main = "Wall to wall")

candidate_metrics <- wall_poly_roads #From previous processing
names(candidate_metrics) <- c("avg", "cov", "std","p10", "p20","p50","p70","p95", "p99","d0","d2","d4","dns")


# Plot candidate metrics
plot(candidate_metrics$p95, zlim = c(0, 30), main = "Candidates")

We calculate PCA values from the lidar metrics located in forested cells using the rasterPCA function from the RStoolbox package.

poly <- st_read("E:/ROMEO/RMF_Pilot_Data/Scripts/RMFinventory/inst/extdata/inventory_polygons/inventory_polygons.shp")
#> Reading layer `inventory_polygons' from data source `E:\ROMEO\RMF_Pilot_Data\Scripts\RMFinventory\inst\extdata\inventory_polygons\inventory_polygons.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 632 features and 111 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 431100 ymin: 5337700 xmax: 438560 ymax: 5343240
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=17 +ellps=GRS80 +units=m +no_defs

#Match CRS of wall-to-wall raster layer
poly <- st_transform(poly, wall_crs@projargs)

#Select FOR polygons and transform to sp
poly_FOR <- poly[poly$POLYTYPE == "FOR", ]
wall_metrics_FOR <- mask(wall_metrics, poly_FOR)

PCA <- RStoolbox::rasterPCA(wall_metrics, nComp = 2, spca = TRUE, maskCheck = FALSE)
# nComp = 2, we return the two first principal components
# spca = TRUE, since metrics have different ranges of values the fucntion will center and scale the metrics
# maskCheck = FALSE, we don't check if some pixels have NA value in one or more layers. If not sure, set to TRUE

# The output of rasterPCA is a list with an element model which contains the PCA model and an element map which contains the map of PCA values

PCA_wall <- PCA$map
plot(PCA_wall)

PCA_model <- PCA$model

We can check the porportion of variance contained in each principal component:

summary(PCA_model)
#> Importance of components:
#>                          Comp.1    Comp.2    Comp.3     Comp.4    Comp.5
#> Standard deviation     2.987499 1.2160491 1.1451346 0.81997529 0.5312166
#> Proportion of Variance 0.686550 0.1137520 0.1008718 0.05171996 0.0217070
#> Cumulative Proportion  0.686550 0.8003019 0.9011737 0.95289369 0.9746007
#>                           Comp.6     Comp.7      Comp.8      Comp.9
#> Standard deviation     0.3919551 0.32657041 0.201570338 0.125814057
#> Proportion of Variance 0.0118176 0.00820371 0.003125431 0.001217629
#> Cumulative Proportion  0.9864183 0.99462201 0.997747439 0.998965068
#>                             Comp.10      Comp.11      Comp.12      Comp.13
#> Standard deviation     0.0843173270 0.0638809609 0.0417904109 2.274839e-02
#> Proportion of Variance 0.0005468778 0.0003139059 0.0001343414 3.980688e-05
#> Cumulative Proportion  0.9995119458 0.9998258517 0.9999601931 1.000000e+00

Use PCA model to get PCA values of an existing set of plots

A network of 182 plots was already established in RMF. The objective of the structural guided sampling was to check if the existing plot network was covering the entire range of structural variability and if not, selecting new plots in underrepresented forest types.

We load the existing set of plots. All the ALS metrics that were calculated to make the PCA model have also been calculated at the plot-level.

Note: Make sure that the name of metrics used to get the PCA model correspond to the name of the plot-level metrics

We can use the predict function to get PCA values of the existing set of plots and off the candidate cells

# Plot-level
plots_metrics <- st_read(system.file("extdata", "plots_metrics.shp", package = "RMFinventory"))
#> Reading layer `plots_metrics' from data source `C:\Users\queinnec.stu\Documents\R\win-library\3.6\RMFinventory\extdata\plots_metrics.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 5 features and 14 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 432109.8 ymin: 5337921 xmax: 437372.2 ymax: 5342820
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=17 +ellps=GRS80 +units=m +no_defs

plots_metrics_coords <- st_coordinates(plots_metrics)
plots_metrics_df <- plots_metrics %>% st_drop_geometry()

plots_PCA <- as.data.frame(predict(PCA_model, plots_metrics_df))[,c(1,2)]
colnames(plots_PCA) <- c("PC1", "PC2")

#Candidate cells

# This is the sricpt that would be ran. Since candidates are just a subset of the wall to wall cells, we could also just mask the PCA_wall layer cells tha are not part of the candidate cells. 
PCA_candidates <- raster::predict(candidate_metrics, model = PCA_model, index = c(1,2),  filename = "path/to/file.tif")

Once PCA values of all forested cells, candidate cells and existing plots have been calculated, it is useful to plot them to visualize their distribution.

df_PCA_wall <- as.data.frame(PCA_wall, na.rm = TRUE, xy = FALSE)
df_PCA_candidates <- as.data.frame(PCA_candidates, na.rm = TRUE, xy = FALSE)

# Get convex hulls 

hulls_wall_idx <- chull(df_PCA_wall$PC1, df_PCA_wall$PC2)
hulls_wall <- dplyr::slice(df_PCA_wall[,c("PC1","PC2")],hulls_wall_idx)

hulls_cand_idx <- chull(df_PCA_candidates$PC1, df_PCA_candidates$PC2)
hulls_cand <- dplyr::slice(df_PCA_candidates[,c("PC1","PC2")],hulls_cand_idx)

ggplot(mapping = aes(x = PC1, y = PC2)) + 
  geom_polygon(data = hulls_wall, colour = "#440154FF", fill ="#440154FF", alpha = 0.2) +
  geom_polygon(data = hulls_cand, colour = "orange", fill = NA) + 
  geom_point(data = df_PCA_candidates, colour = "orange") + 
  geom_point(data = plots_PCA) + 
  theme_bw() +
  coord_equal() +
  theme(panel.grid = element_blank(),
        legend.position = "right")

Stratification

Once PCA values have been calculated we can stratify the PC1 / PC2 feature space. The approach chosen for the RMF was to stratify using equal intervals.

From the previous plot, we can decide to stratify the feature space into 7 equal interval on the PC1 axis and 5 intervals on the PC2 axis (7 x 5 grid). We can use the function getPCAstrata to automatically create the stratification:

# Determine the number of breaks for the first and second features (PC1 and PC2) 
strata <- RMFinventory::getPCAstrata(PCA_layer = PCA_candidates,
                                  nbreaks = c(8, 6), #Since we want a 7 x 5 grid we need 8 x 6 breaks
                                  summary = TRUE)


#getPCAstrata returns a list of three objects

strata_candidates <- strata$strata_layer
strata_candidates
#> class      : RasterLayer 
#> dimensions : 277, 373, 103321  (nrow, ncol, ncell)
#> resolution : 20, 20  (x, y)
#> extent     : 431100, 438560, 5337700, 5343240  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
#> source     : memory
#> names      : strata 
#> values     : 12, 74  (min, max)
plot(strata_candidates)

breaks <- strata$breaks
breaks
#> $PC1
#> [1] -5.93029133 -3.92324124 -1.91619115  0.09085893  2.09790902  4.10495911
#> [7]  6.11200920  8.11905929
#> 
#> $PC2
#> [1] -5.10888294 -3.40139913 -1.69391531  0.01356851  1.72105233  3.42853614

strata$summary

strata$matrix
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,]   11   21   31   41   51   61   71
#> [2,]   12   22   32   42   52   62   72
#> [3,]   13   23   33   43   53   63   73
#> [4,]   14   24   34   44   54   64   74
#> [5,]   15   25   35   45   55   65   75

We can use the breaks returned by the getPCAstrata function to get the strata of the existing set of plots:

strata <- getPCAstrata(PCA_layer = plots_PCA,
                             breaks = breaks, # From the precedent call to getPCAstrata
                             summary = TRUE)
#> breaks have been provided and won't be recalculated

strata_plots <- strata$strata_layer
strata_plots
#> [1] 44 24 52 54 51
#> 35 Levels: 11 < 12 < 13 < 14 < 15 < 21 < 22 < 23 < 24 < 25 < 31 < ... < 75

strata$summary

Make a new plots with break lines between strata

ggplot(mapping = aes(x = PC1, y = PC2)) + 
  geom_polygon(data = hulls_wall, colour = "#440154FF", fill ="#440154FF", alpha = 0.2) +
  geom_polygon(data = hulls_cand, colour = "orange", fill = NA) + 
  geom_point(data = df_PCA_candidates, colour = "orange") + 
  geom_point(data = plots_PCA) + 
  geom_vline(xintercept = breaks$PC1, linetype = "dashed") + 
  geom_hline(yintercept = breaks$PC2, linetype = "dashed") + 
  theme_bw() +
  coord_equal() +
  theme(panel.grid = element_blank(),
        legend.position = "right")

Select new plots

The function sampleCells that performs the sampling requires the number of cells to sample for each strata. This can be determined based on the previous feature space plots, total number of plots than can be selected, number of exisiting plots, how many of them should be re-measured etc.

The following figure illustrates the sampling process:

Sampling process diagram

toSample <- read.csv(system.file("extdata","cells_to_sample.csv", package = "RMFinventory"), stringsAsFactors = F)
toSample

existing_plots <- data.frame(plotID = plots_metrics$plotID, 
                             x = plots_metrics_coords[,1], 
                             y = plots_metrics_coords[,2], 
                             strata = strata_plots)
existing_plots

new_plots <- RMFinventory ::sampleCells(strata_layer = strata_candidates,
                         matrix_strata = strata$matrix, 
                         existing_sample = existing_plots, # You can provide a set of exisitng plots or output of previous call to sampleCells and these cells won't be sampled again.
                         toSample = toSample, 
                         mindist = 150, 
                         message = T)
#> Strata 13: 1 cells to sample
#> Strata 14: 2 cells to sample
#> Strata 22: 1 cells to sample
#> Strata 23: 3 cells to sample
#> Strata 24: 3 cells to sample
#> Selecting in extended clusters ...
#> Strata 25: 0 cells to sample
#> Strata 31: 1 cells to sample
#> Selecting in extended clusters ...
#> Selecting isolated cells ...
#> Strata 32: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 33: 5 cells to sample
#> Strata 34: 3 cells to sample
#> Selecting in extended clusters ...
#> Strata 35: 0 cells to sample
#> Strata 41: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 42: 2 cells to sample
#> Selecting in extended clusters ...
#> Strata 43: 5 cells to sample
#> Strata 44: 3 cells to sample
#> Strata 45: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 51: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 52: 2 cells to sample
#> Selecting in extended clusters ...
#> Strata 53: 4 cells to sample
#> Selecting in extended clusters ...
#> Strata 54: 5 cells to sample
#> Strata 55: 2 cells to sample
#> Strata 61: 0 cells to sample
#> Strata 62: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 63: 2 cells to sample
#> Selecting in extended clusters ...
#> Strata 64: 2 cells to sample
#> Selecting in extended clusters ...
#> Strata 65: 2 cells to sample
#> Selecting in extended clusters ...
#> Strata 72: 0 cells to sample
#> Strata 73: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 74: 1 cells to sample
#> Selecting in extended clusters ...
#> Strata 75: 0 cells to sample

# There might be a warning like this: 
#no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
# If think safe to ignore but need to look more into that
new_plots 

new_plots <- st_as_sf(new_plots, coords = c("x", "y"))

plot(strata_candidates)
plot(st_geometry(new_plots), add = T)

new_plots_metrics <- extract(wall_metrics, new_plots)
new_plots_PCA <- as.data.frame(predict(PCA_model, new_plots_metrics))[,c(1,2)]
colnames(new_plots_PCA) <- c("PC1", "PC2")


ggplot(mapping = aes(x = PC1, y = PC2)) + 
  geom_polygon(data = hulls_wall, colour = "#440154FF", fill ="#440154FF", alpha = 0.2) +
  geom_polygon(data = hulls_cand, colour = "orange", fill = NA) + 
  geom_point(data = df_PCA_candidates, colour = "orange") + 
  geom_point(data = new_plots_PCA) + 
  geom_vline(xintercept = breaks$PC1, linetype = "dashed") + 
  geom_hline(yintercept = breaks$PC2, linetype = "dashed") + 
  theme_bw() +
  coord_equal() +
  theme(panel.grid = element_blank(),
        legend.position = "right")

Field inventory data

Forest attributes using allometric equations

To be completed

Calculate plot-level ALS metrics

Plot-level ALS metrics need to be calculated in order to build predictive models of forest attributes. These task can be performed in various softwares/packages. Here, we present two of them: LAStools and the lidR R package.

The lidR package offers more flexibility to calculate user-defined metrics. At the plot-level, computing time difference between lidR and LAStools is not significant. However, ALS metrics also need to be calculated wall-to-wall in order to generate forest attribute maps once predictive models have been developed on plot data. Over larger areas, LAStools will generally be more efficient than lidR.

In the end, the choice between lidR and LAStools is a trade-off between flexibility in ALS metrics and computational efficiency when calculating wall-to-wall maps.

Using the lidR package

The functionalities of the lidR package are well described in the following wiki: https://github.com/Jean-Romain/lidR/wiki.

The Area-based approach from A to Z is a great resource to have a look at for ABA.

Clip to plots extent

lastiles <- lidR::readLAScatalog("path/to/ALS tiles")

plots <- st_read("path/to/plots_coordinates.shp")

opt_output_files(lastiles) <- "path/to/output/plot_{PlotID}" # Here plotID is a field of plots_coordinates.shp

plot_radius = 11.28
plots_las <- lasclip(lastiles, plots, radius = plot_radius)

Calculate plot_level metrics

To be completed

Using lastools

To be completed

Clip to plots extent

To be completed

Calculate plot-level metrics

It is possible to execute operative system commands through R using the system function by passing it a character containing the command.

files <- "path/to/lasfiles"
odir <- "path/to/output/file.csv"
lastools_cmd <- paste0("lascanopy -i ", files, "*.las -files_are_plots -height_cutoff 1.3 -cover_cutoff 2 -p 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 99 -avg -cov -dns -std -kur -ske -d 2 5 10 15 20 25 -names -o ", odir)
system(lastools_cmd)

Make predictive models

The makeModels function implemented in the package builds random forest predictive models for each forest attribute and assess their accuracy using a k-fold cross-validation.

# Open plot-level ALS metrics

als_metrics <- read.csv(system.file("extdata","plot_metrics.csv", package = "RMFinventory"))
head(als_metrics)

toKeep <- colnames(als_metrics)[colnames(als_metrics) %in% names(wall_metrics)]

als_metrics <- als_metrics[,c("plotID", toKeep)]

pred_names <- colnames(als_metrics)[!colnames(als_metrics) %in% "plotID"]
pred_names
#>  [1] "avg" "std" "p10" "p20" "p50" "p70" "p95" "p99" "d0"  "d2"  "d4" 
#> [12] "cov" "dns"

# Open plot-level forest attributes 

plot_data <- read.csv(system.file("extdata","tblPlotLevelSummary.csv", package = "RMFinventory"))
head(plot_data)

# Combine field data and ALS metrics
combine_dat <- merge(als_metrics, plot_data, by = "plotID")

# Attributes to model
attNames <- c("top_height", "ba_ha", "V_ha") 
# Title to add to plots (optional)
titles <- list(ba_ha =  expression("Basal area (m"^"2"*"/ha)"), 
               top_height = "Top height (m)", 
               V_ha = expression("Whole stem volume (m"^"3"*"/ha)"))

# Run modeling: k-fold cross validation, random forest regression models and accuracy assessment
rfList <- RMFinventory::makeModels(dat = combine_dat,
                                   attNames = attNames,
                                   preds = pred_names, 
                                   k = 5, 
                                   titles = titles,
                                   saveModel = F,
                                   saveFigure = F)
#> Loading required package: lattice
rfList$accuracy$`all folds`

rfList$accuracy$summary

# The accuracy measures reported on the plot are currently not extracted from rfList$accuracy$summary as they should be. They are calculated based on all plotted points. 
# The mean (and std) of accuracy measures of the cross-validation are better indicators of the "true" accuracy assessment. 
rfList$top_height$plot

rfList$ba_ha$plot

rfList$V_ha$plot

Calculate wall-to-wall ALS metrics

# Wall to wall metrics (already opened at the beginnign of the vignette)
wall_metrics <- brick(system.file("extdata", "wall_metrics_small.tif", package = "RMFinventory"))
names(wall_metrics) <- c("avg", "cov", "std","p10", "p20","p50","p70","p95", "p99","d0","d2","d4","dns")

Get maps of forest attributes

top_height_wall <- raster::predict(wall_metrics, rfList$top_height$finalModel)
ba_ha_wall <- raster::predict(wall_metrics, rfList$ba_ha$finalModel)
V_ha_wall <- raster::predict(wall_metrics, rfList$V_ha$finalModel)

plot(top_height_wall, main = titles$top_height)

plot(ba_ha_wall, main = titles$ba_ha, zlim = c(0, 70))

plot(V_ha_wall, main = titles$V_ha, zlim = c(0, 500))